/*
** LuaFileSystem
** Copyright Kepler Project 2003 - 2017 (http://keplerproject.github.io/luafilesystem)
*/

#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
#include <Eigen/SparseLU>
#include <log/log.h>
#include <Spectra/GenEigsSolver.h>
#include <Spectra/SymEigsSolver.h>

#include <Spectra/SymGEigsShiftSolver.h>
#include <Spectra/SymGEigsSolver.h>
#include <Spectra/MatOp/SparseGenMatProd.h>
#include <Spectra/MatOp/SparseSymMatProd.h>
#include <Spectra/MatOp/SparseCholesky.h>
#include <Spectra/MatOp/SparseRegularInverse.h>

#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <Eigen/Eigenvalues>
#include <Spectra/SymGEigsSolver.h>
#include <Spectra/MatOp/DenseSymMatProd.h>
#include <Spectra/MatOp/SparseCholesky.h>

//------------------------------------------------------------------------------//
//要求： A是实对称矩阵,下三角矩阵，三角数据被忽略
//   ： B是实对称矩阵,下三角矩阵，三角数据被忽略
template<class T, class T2>
int solve_gsheig_spectra(Eigen::SparseMatrix<T> &A, Eigen::SparseMatrix<T> &B,
		Eigen::Matrix<T2, Eigen::Dynamic, 1> &evalues,
		Eigen::Matrix<T2, Eigen::Dynamic, Eigen::Dynamic> &evecs,
		int nbrEigenvalues) {
	using namespace Spectra;

	SparseSymMatProd<double> op(A);
	SparseCholesky<double> Bop(B);
	SymGEigsSolver<SparseSymMatProd<double>,SparseCholesky<double>,GEigsMode::Cholesky>
				geigs(op, Bop, nbrEigenvalues, nbrEigenvalues+3);
	geigs.init();
	int nconv = geigs.compute(SortRule::SmallestMagn);
	if (geigs.info() == CompInfo::Successful){
		evalues = geigs.eigenvalues();
		evecs = geigs.eigenvectors();
	}else{
		return -1;
	}
	return 0;
}
//------------------------------------------------------------------------------//

